Method and apparatus for reducing the effects of motion in an image

ABSTRACT

A method and apparatus for reducing the effects of motion of a subject during image acquisition is presented. The method is particularly suitable for correcting for motion during magnetic resonance imaging. The method involves comparing an image with a reference image in the spatial frequency domain to determine the extent of any rotational and/or translational motion and then apply an appropriate correction. Any motion is modelled by comparison with the reference image but the corrected image is produced using information from the corrupted image alone.

[0001] This invention relates to the field of image processing and inparticular to image processing techniques whereby artefacts introducedinto an image during acquisition, especially through motion of theobject being imaged, are corrected for in the generation of an outputimage.

[0002] Magnetic Resonance Imaging (MRI) is a widely used technique formedical diagnostic imaging. In a conventional MRI scanner, a patient isplaced in an intense static magnetic field which result in the alignmentof the magnetic moments of nuclei with non zero spin quantum numberseither parallel or anti-parallel to the field direction. Boltzmanndistribution of moments between the two orientations results in a netmagnetisation along the field direction. This magnetisation may bemanipulated by applying a radio frequency (RF) magnetic field at afrequency determined by the nuclear species under study and the strengthof the applied field. In almost all cases for medical imaging thespecies studied is the nucleus of the hydrogen atom present in the bodymainly in water molecules, and the RF pulses are applied at the resonantfrequency of these water protons. However it is also possible to imageother nuclear species and hydrogen bound in other species.

[0003] The energy absorbed by nuclei from the RF field is subsequentlyre-emitted and may be detected as an oscillating electrical voltage, orfree induction decay signal, in an appropriately tuned antenna. Morecommonly a further RF pulse or magnetic field gradient is used topostpone signal acquisition and generate a spin-echo or gradientecho-echo signal.

[0004] Spatial information is encoded into the echo signal by virtue ofadditional linearly varying magnetic fields, known as magnetic fieldgradients, applied during or prior to echo acquisition. The principle ofspatial encoding is that in the presence of the field gradient the netfield experienced by a given nuclear moment, and hence it resonantfrequency, is a function of position with the scanner. When a gradientis applied during echo acquisition the received signal contains a rangeof frequency components representing nuclei at different locations alongthe gradient direction. Fourier transformation of this signal yields aone-dimensional projection through the patient. This technique is knownas frequency encoding. Two dimensional encoding requires use of anadditional gradient, applied perpendicular to the frequency encodingaxis, known as the phase encoding gradient. This gradient is applied fora short time prior to data acquisition. The acquisition process isrepeated a number of times, perhaps 256 or 512 times, using phaseencoding gradients of different strengths. Simultaneous frequency andphase encoding yields a two-dimensional data set which, when subjectedto a two-dimensional Fourier transform, provides the required image.This array of data exists in what is known as k-space and is the Fouriertransform of the image space. The effect of the phase encoding gradientis to move the start of the data acquisition to a particular locationalong one axis in k-space (dependent on the gradient strength) whilstfrequency encoding represents a sweep through k-space parallel to theother axis. Each of these sweeps in known as a “shot” or “view”.

[0005] Spatial localisation in the third dimension may be achieved usingan additional phase encoding gradient, or more commonly by using agradient and narrowband RF pulse to restrict the initial perturbation ofnuclear moments to a single tomographic slice. This principle can bereadily extended to multislice MRI.

[0006] In conventional MRI a single phase encoding view is acquiredafter each RF excitation. However, faster imaging sequences now exist inwhich further RF pulses and phase encoding gradients are used to acquirea train of differently encoded echoes after each excitation. Theseechoes traverse several lines of k-space and reduce scanning time by afactor equal to the echo train length. In the extreme case single shotecho planar imaging (EPI) techniques cover the whole of two dimensionalk-space in a single acquisition lasting less than 100 ms, althoughspatial resolution and image quality may be significantly compromised.

[0007] Patient movement during the acquisition of MRI images results indegradation of the images that can obscure the clinically relevantinformation. Each readout period takes a few milliseconds whereas thetime interval between readouts might be between 100 and 4000 ms. Themajority of blurring and ghosting artefacts caused by patient motion aredue to motion between lines in k-space, rather than motion during asingle readout.

[0008] Movement leads to phase errors between lines of k-space which, inthe resulting image, appear as blurring and hosting along the phaseencode direction. These phase errors can result from translationalmovement and rotational movement. Translations of the patient in thereadout direction result in a frequency dependent phase shift in eachline of k-space. Rotations in the spatial domain are also rotations ink-space and result in k-space changes that are more complicated functionof position in k-space.

[0009] Various techniques have been employed to correct for imageartefacts introduced into an image through motion. However most of thetechniques known for correcting for patient motion involve a modifiedsignal acquisition technique which may involve additional scans or evenadditional equipment.

[0010] International Patent Application WO98/01828 discloses a techniquefor reducing the effect of motion induced artefacts in an image usingpurely post data gathering signal processing effects. In the techniquedescribed therein the data is manipulated to counteract possiblemovement induced artefacts and the manipulated data compared using afocus condition to see if the image quality is improved. This techniquecan involve a large amount of processing however. Furthermore the methodmay involve the grouping of k-space lines to more accurately determinemotion parameters, however this grouping can decrease the temporalresolution of the motion found.

[0011] Another method of correcting for motion induced image artefactsin the method of Projection onto Complex Sets (POCS) Hedley M, Hong Yand Rosenfeld D. “Motion Artifact Correction in MRI using generalizedprojections” IEEE Trans. Med. Imag., 10:40-46, 1991. This is a methodwhereby a good quality image is used to form a binary mask. The maskdefines the tissue-air boundary, i.e. outside the mask there should beno signal. Motion induced artefacts in the acquired image cause apparentsignal in the air. The POCS method sets all outside the mask in theacquired image to black. The image data is then Fourier transformed tok-space. A new complex k-space is formed from the modulus of themeasured data and the phase of the estimation from the previous step.This new k-space is Fourier transformed to the image domain and theprocess iterates. This method however involves a large amount of Fouriertransformation as the process iterates and hence involves a large amountof computational effort and hence time. Further the method requires thespatial alignment of the binary mask with the acquired image beforeprocessing which is not always possible to achieve.

[0012] A different method of correcting for motion induced artefacts indiffision-weighted MRI is discussed in ‘Correction of Motional Artifactsin Diffusion-Weighted Images Using a Reference Phase Map’, Ulug et al.Magnetic Resonance in Medicine, Vol. 34, no. 2, 1995, pp476-480. Themethod described therein takes a measurement with no diffusionweighting, or a very low weighting, and uses this information toconstruct a phase map. The phases from this phase map are then used toreplace the phases of higher weighted data sets to eliminate a motioninduced phase component and form an adjusted data set. Images can thenbe formed from this motion data set.

[0013] The method described by Ulug et al. concerns diffusion-weightedMRI scans, which are very susceptible to motion, and is not applicableto all MRI images. Further the method described would not be capable ofcorrecting for rotations of the object during imaging.

[0014] The present invention seeks to provide an alternative method forcorrecting for motion induced artefacts in an image.

[0015] Thus according to the present invention there is provided

[0016] A method of producing an image of a scanned object corrected forartefacts introduced by unwanted motion of said object during the scancomprising the steps of:

[0017] taking a scan k-space data set derived from the object scan,

[0018] taking a reference k-space data set of the object,

[0019] comparing the scan k-space data set and reference k-space dataset to determine the extent of any motion during the scan in the scank-space data set,

[0020] using the scan k-space data set and information on the extent ofmotion of the scan k-space data set to produce a motion corrected scank-space data set, and

[0021] using the image corrected scan k-space data set to produce animage.

[0022] The reference k-space data set is a motion free (or as near tomotion free as possible in the circumstances) k-space data setcorresponding to an image of the object. By taking a motion free dataset the motion leading to the corrupted image data set can be modelled.The data from the scan k-space data set alone is then used to create animage corrected on the basis of the modelled motion. As the two datasets are compared in the k-space or spatial frequency domain there is noneed for repeated Fourier transforms of the whole data set. As themethod corrects the data of the scan k-space data set for determinedmotion rather than replace it with other data the information from theactual scan is preserved.

[0023] The method preferably comprises the step of dividing the scank-space and reference k-space data sets into sections in k-space and thestep of comparing the scan k-space data set and reference data setcompares the data sets section by section to determine the motion. For2D imaging each section could conveniently be a line of phase encodedata taken by the scanner. For 3D imaging each section is a plane ofk-space data.

[0024] Conveniently the step of comparing the two data sets comprisesthe step of comparing the data sets to determine the extent of anyrotation and also the step or comparing the data sets the extent of anytranslation. Motion that causes an expansion or contraction in the imagemight also be corrected by such a method.

[0025] Preferably the step of comparing the data sets includes the stepof comparing selected points in the data sets to determine the extent ofany rotational motion. Conveniently the k-space moduli of selectedpoints can be used for the comparison. The modulus of the k-space datais independent of phase and as such is not affected by translationalmotion. Rotations of the object being imaged cause rotations in k-spacewhereas translations impart phase ramps to the k-space data. Theselection of the points may involve separating the k-space data intoblocks. Certain points within a block could be selected, for instance byselecting the points with the largest modulus value in each block. Theseparation into blocks ensures that a reasonable coverage of k-space isachieved to give good angular resolution.

[0026] The extent of rotational motion can be determined by iterativelyconsidering possible angular motions of the selected points in thereference k-space data and comparing the transformed reference data setwith the scan k-space data set using appropriate focus criteria todetermine the extent of rotational motion. The step of transforming thereference k-space data set is preferably carried out using linearinterpolation although other interpolations such as cubic or splinecould be used.

[0027] Conveniently the step of comparing the data sets to determine theextent of translational motion are determined using a method based on across spectrum method such as described by Reddy B S and Chatterji E N,“An FFT-Based Technique for Translation, Rotation and Scale-InvariantImage Registration” IEEE Trans. Imag. Proc., 5:1266-1271, 1996. Inperforming the cross spectrum analysis k-space may be zero padded or maybe padded with complex values having a modulus equal to one and randomphases in order to improve the resolution. The cross-spectrum methoddetects constant and linearly varying phase differences between the twodata sets. The two data sets are combined to form a function K_(c) aswill be described below. The Fourier Transform of K_(c) gives a peak.The shift of the peak from centre is related to the phase slope andgives an indication of the amount of motion in the FE direction. Thephase of the peak reveals any phase difference between the two data setsand reveals the extent of motion in the PE direction.

[0028] After determining the extent of translational motion the scank-space data set may be corrected for the determined motion. The phaseof this corrected data set could then be compared to the phase ofreference k-space data to determine sub-pixel motion in the frequencyencode direction.

[0029] Where any point is known to be unreliable in either the referencek-space data set or scan k-space data set the corresponding point inK_(c) will be unreliable. Such a point could be replaced with a pointhaving a modulus equal to one and random phase to prevent the unreliablepoint affecting the determination of motion.

[0030] The cross-spectrum based method may include a measure ofconfidence. In the ideal case the result of the Fourier Transform ofK_(c) would be a single peak. In reality differences between thereference and scan k-space data sets due to noise, different acquisitioncriteria etc. will have an effect. This will result in the FT of K_(c)being a peak with other values that are non-zero. The ratio of theheight of the peak to the average of the other values could be used as ameasure of confidence, the higher the ratio the more reliable thedetermination of motion.

[0031] The steps of taking the reference k-space data set and scank-space data set may include the step of aligning both the scan k-spacedata set and reference k-space data set so that the maximum (DC) pointcoincides with the centre of the data grid. Further the method caninclude the step of globally scaling the two data sets so that at themaximum (DC) point the two data sets have the same complex value.

[0032] Note that the k-space point with the maximum modulus value isassumed to be the point that corresponds to zero spatial frequency. Zerospatial frequency is sometimes called DC.

[0033] Further the method may include the step of filtering both theacquired image and reference image to introduce a certain amount ofblurring. Applying a filter to the images blurs the correspondingk-spaces and can aid the determination of rotational motion using linearinterpolation. Blurring the k-space means that it varies less rapidlyand so interpolation between points is easier.

[0034] The scan k-space data set may be a data set obtained from amagnetic resonance imaging (MRI) scanning means.

[0035] The reference k-space data set can be obtained in any number ofways. Many scanning protocols acquire two or more data sets and averagein order to provide better to signal to noise ratio. In such cases oneof the data sets could be taken as the reference image to correct forany motion in the other(s). The best data set could be determined byusing a measure of the entropy of the image, i.e. the degree of disorderin the image.

[0036] Alternatively, where the object is being imaged successively inconjunction with introduction of a contrast agent one of the imagescould be used as the reference k-space data set.

[0037] Alternatively the reference k-space data set could be taken froma previously acquired scan. In monitoring of a progressive disorder saythere could be a need for repeated scans on the same subjects atdifferent times. As the present invention does not use any k-spaceinformation from the reference data set but only determines the extentof any motion for the scan k-space data set a previous image can be usedand the technique should be reasonably robust against differencesbetween images such as differing contrast, noise or even features suchas an enhancing tumour. Where no reference data is available a referenceimage could be fabricated from the available data.

[0038] Rather than correct a whole image, regions of the image can bemasked permitting correction on a localised region only. In this casethe method includes the steps of selecting a desired part of the imageand blanking the non required areas of the image before Fouriertransforming the image to produce the scan k-space data set.

[0039] In another aspect of the invention there is provided an imagingsystem comprising a scanning means for scanning an object and driving ascan k-space data set from signals received from the object and aprocessing means for manipulating the scan k-space data set to producean image of the scanned object corrected for artefacts introduced byunwanted motion of said object during the scan wherein the processingmeans is adapted to take the scan k-space data set derived from theobject scan, take a reference k-space data set of the object, comparethe scan k-space data set and reference k-space data set to determinethe extent of any motion in the scan data set, use the information onthe extent of motion of the scan k-space data set to produce a motioncorrected scan k-space data set, and use the motion corrected scank-space data set to produce an image.

[0040] The invention will now be described by way of example only withreference to the following drawings of which;

[0041]FIG. 1 shows a schematic diagram of a magnetic resonance imagingsystem;

[0042]FIG. 2 shows a functional diagram of the operation of the FIG. 1system;

[0043]FIG. 3 shows a flow diagram of an image correction routineperformed according to the present invention;

[0044]FIG. 4 shows an MRI image, a reference image and a motioncorrected image using the method of the present invention.

[0045]FIG. 5 shows orthogonal MRI images for an acquired image, areference image and a motion corrected image using the presentinvention.

[0046] Referring to FIG. 1 there is shown a schematic diagram of amagnetic resonance imaging system 10. The system 10 incorporates amagnetic resonance imaging scanner 12 of conventional type. The scanner12 has a superconducting or resistive main magnet 20 which generates amagnetic field sufficiently strong to cause a net alignment along thefield direction of atomic nuclei within a patient. The scanner 12 alsoincludes shim coils 22 in order to correct for undesired inhomogeneitiesin the magnetic field of the main magnet 20. The magnetic field producedby the shim coils 22 is controlled by a shim coil power supply unit 24.

[0047] The resonance frequency of particular atomic nuclei ischaracteristic of the nucleus and the strength of the applied magneticfield. In order to provide spatial information a magnetic field gradientis generated by gradient coils such as coils 26. Gradient coils areoften arranged to generate gradient fields in three orthogonaldirections. The magnetic fields generated by the gradient coils arecontrolled by a gradient coil power supply unit 28. In order to generatea signal from the atomic nuclei or the patient a radio-frequencymagnetic pulse is generated by transmit coil 30. This pulse ‘flips’ theangle of the nuclear spins within a certain patient slice of volume.These excited spins or magnetisations then induce a current in thereceive coil which may be the same coil as the transmit coil 30. Thecoil 30 is connected to a transmit unit 32 and a receive unit 34, eachof which also receives signals from a frequency source 36.

[0048] The system 10 includes a controlling computer 38 which controlsthe operation of the components of the system 10. The computer 38controls the gradient coil power supply unit 28 in the form of gradienttiming, magnetic field strength and orientation control. In addition,the computer receives signals from the receive unit 34 together withtransmitter timings.

[0049] In order to form an image of the organs of a patient, the patientis inserted into the system 10 and a series of measurements are takenwith different combinations of static and/or varying gradient fields.The signals from the tissue of the patient depend on the tissue'sproperties, the magnetic field gradient strengths, gradient orientationsand timings with respect to the applied radio frequency pulses. Thevarying gradients code the received signal's phase, frequency andintensity. The received signals as a function of time form an orderedset which is stored in memory in the computer 38 for subsequentprocessing.

[0050] In a subsequent signal processing stage a Fourier transform maybe performed on the ordered set of received signals, with the modulus ofthe transform being used to assign the signals to a grey scale in orderto form an image. The set of received signals is said to exist ink-space.

[0051] In a conventional MRI if a patient moves during the acquisitionof data the received signal is affected and part of the k-space signalis corrupted. Because of the way the image is reconstructed this motionaffects the whole image, causing blurring and/or ghosting artefacts inthe final image.

[0052] Referring now to FIG. 2 there is shown a functional block diagramof the operation of the system 10. The computer 38 controls and receivesinformation from the scanner 12 and uses this information to generate animage on display 50. This image is an initial reconstructed image. If anoperator of the system 10 considers that the initial image is corruptedan additional signal processing routine is selected. Alternatively thefurther signal processing could occur automatically. In either case, thestored image data is processed to reduce the effects of the patient'smotion.

[0053] The method of processing the image is as follows with referenceto FIG. 3 which shows a flow chart 100 of the stages in an imageprocessing routine according to the present invention.

[0054] The acquisition time is divided into nodes. For 2D images a nodecorresponds to a phase encode line. For 3D images a node is a plane ofk-space data.

[0055] The axes of the system are defined as follows: X is the readout,or frequency encode direction. For 2D imaging Y is the phase encodegradient direction and Z is the through slice direction. In 3D imaging Yis the slow phase gradient encode direction and Z is the fast phaseencode gradient direction. Therefore data points are acquired with k_(x)varying most rapidly, then k_(z) and most slowly in k_(y). In 3D imaginga node is a plane of constant k_(y).

[0056] The image processing is carried out on image data 102, which maybe obtained, for instance, as sequential k-space lines although otheracquisition techniques are possible. A reference image data set 104 isalso used.

[0057] Initially a degree of pre-processing of the data sets may beperformed to aid in the comparison. In this case the first stage may bean alignment stage 106 where both k-spaces are aligned such that themaximum (DC) value k-space data point for each set is at the centrepixel for both data sets.

[0058] The two data sets may then be scaled in a normalisation stage 108to ensure that the complex value at DC is the same in each data set.This is preferably done on the reference data set 104.

[0059] A filtering stage 110 is then employed to filter the actual imageand the reference image. This filter introduces a degree of blurringinto the k-space which aids in the comparison stages later. A suitablefiltering function is to multiply the pixels by a Gaussian with value 1at the image centre falling away towards the image edge. One filter thathas been used is a Gaussian with a variance equal to (0.3125 n)² where nis the number of points along that k-space dimension.

[0060] By filtering in such a manner k-space is blurred (and the imageis bright in the middle and gets darker towards the edges). This aidsthe linear interpolation step later.

[0061] If required a masking stage 112 can be applied where pixels notof interest can be set to zero.

[0062] The filtered and suitably masked images then undergo a Fouriertransform stage 114 to produce the k-space image data set and k-spacereference data set. The filtering and masking is performed on the actualimage hence the need to Fourier transform into k-space.

[0063] In order to ensure that appropriate angular resolution isachieved k-space is then partitioned into various blocks in apartitioning stage 116. Splitting k-space into blocks can enable anincrease in speed of the algorithm but can result in a loss of accuracyand this partitioning can be omitted by making every k-space point aseparate block. For a 2D image there may be a number of blocks spacedevenly along the line. If there were 256 points in a line choosing 32blocks would result in 8 points per block although other combinationscould be used. In 3D the plane is divided into a certain number ofblocks in each dimension of the plane.

[0064] Next, at the node selection stage 118 each node is selected inturn. For each node the coordinates, r₁ of the points with the largestmodulus value in each block are determined in stage 120 although otherselection criteria could be utilized.

[0065] An angular transformation T_(a) corresponding to a rotation ofinterest is selected, different transformation being selected in turn ina selection stage 122. The modulus of the positions in the referencek-space K_(r) corresponding to the coordinates of the selected pointstransformed by the angular transformation T_(a) are then determined. Inthis stage 124 linear interpolation techniques, as will be wellunderstood to one skilled in the art, are sufficient when the filteringstage 110 has been applied.

[0066] The angular transformation T_(a) that most accurately representsthe rotational motion can be determined by finding which angulartransformation minimises a suitable focus criteria.

[0067] A suitable focus criteria is;$f_{c} = {\sum\limits_{i}\frac{\left( \left. {{K_{r}\left( {T_{a}\left( r_{i} \right)} \right.} - {{K_{m}\left( r_{i} \right)}}} \right)^{2} \right.}{{{K_{r}\left( {T_{a}\left( r_{i} \right)} \right)}}^{2} - {{K_{m}\left( r_{i} \right)}}^{2}}}$

[0068] where K_(m)(r_(i)) is the unmodified image data for the selectedpositions. When T_(a)(r_(i)) falls outside the measured referencek-space it should be excluded from the focus criteria sum.

[0069] This focus criteria is suitable for nodes close to the central(C) node. For nodes further away say more than 60 nodes from the central(DC) node the focus criteria sum can be changed so that the numeratorsand denominators are summed separately before taking the ratio. Thenormalisation of the focus criterion can be changed to alter therelative contribution from k-space points with different modulus values.Near the centre where the moduli are larger than the edges each pointcarries the same weight which gives the best angular resolution. For theouter lines however, which have lower moduli, noise is more of a problemand changing the normalisation permits a trade off between noise andangular resolution.

[0070] As an alternative, instead of interpolating points within onedata set the relevant line from a previously rotated and stored data setmight be used in the evaluation of the focus criterion.

[0071] Once all the angular transformations have been checked and theminimum focus criteria determined for that node the next node can betried.

[0072] Once all nodes have been processed a rotation correcting stage isundertaken where, for instance, a shearing algorithm is applied, as iswell known, to apply the rotations found to the k-space data set. Othermethods for correcting such as interpolation could also be appliedhowever.

[0073] Next the translations are corrected starting with again with anode choice step 130. In an alternative embodiment however thetranslational correction for a particular node could be determinedimmediately after a rotational correction for a node has been performed.It is not necessary to wait for all nodes to undergo rotationalcorrection

[0074] The translations are determined using a cross spectrum basedmethod. If K_(r) represents all the k-space points in a line (or a planein 3D) of the reference data at a given node and K_(m) the same in theimage data then the position of the peak in the Fourier transform ofK_(c) gives the translation in the defined x-direction (and z in 3D).K_(c) is given by the equation;$K_{c} = \frac{K_{r}^{*}K_{m}}{{K_{r}{K_{m}}}}$

[0075] The shift of the peak of the Fourier transform of K_(c) from thecentre of the array gives the x and z directions (just x in 2D) and they-direction can be determined from the phase p of the peak. The ydisplacement Y is given by; $Y = \frac{p}{2\pi \quad l}$

[0076] where 1 is the number of the phase encode line, being 0 at thecentre (DC) phase encode line.

[0077] To improve the resolution by which x and z displacements arefound, K_(c) may be zero padded. Alternatively, K_(c) may be padded withcomplex values that have modulus one and random phases. The latter mayenable the peak in the Fourier transform of K_(c) to be more clearlydistinguished and accurately located.

[0078] Further if any points in either K_(r) or K_(m) are known to beunreliable, the corresponding point in K_(c) will also be unreliable. Insuch a case the unreliable point could be replaced with a point havingmodulus one and random phase. A modulus of one is chosen as every pointof K_(c) will have a modulus of one due to the inherent normalisationpresent in K_(c). A point may be determined as unreliable by some otherprocessing step such as data averaging or calculation of k-space frommultiple receiver coils.

[0079] The translation information can then be used in a correctionstage 134 to correct the data for that node before moving onto the nextnode.

[0080] Where the cross-spectrum based method as described above is usedwithout data padding the method determines the translational motion inthe frequency encode direction correct to the nearest pixel. In analternative embodiment (not shown) the method can be run as describedabove using the reference data set and the image data set. The imagedata set can then be corrected for the determined motion and thisintermediate corrected data then compared with the reference data set.The remaining phase ramp can then be determined to give sub-pixel shiftsin the fe direction. The phase ramp may be determined using a robustline fitting algorithm on the phase of the ratio of the K_(r) tocorrected K_(m).

[0081] In this way the acquired image may be compared to the referenceimage and any rotational and/or translational movement of the acquiredimage can be corrected. Carrying out the rotational comparison on thek-space modulus data has the advantage of being invariant to anytranslational motion. Also the centre of rotation in the image domaindoes not need to be known as k-space always rotates about its centre.

[0082] A consequence of this method is that the reference and correctedimage are spatially aligned. Further no prior spatial alignment of thereference and data image sets is required unlike the method Projectiononto Convex Sets.

[0083] Also although all the information in the reference image is usedto inform the motion correction required for the image data all of theactual data points used in the corrected image come from the image data.

[0084] The reference image could be a previously acquired scan. Somepatients with progressive disorders have several scans and a historicalscan could be used. As none of the actual reference image data is usedin the corrected image the method can be used with slight differencesbetween the images such as might be encountered between scans atdifferent stages of a slowly progressing disorder.

[0085] Alternatively where multiple images are taken, for instance withincreasing contrast agent a particular image could be taken as thereference image. Further some imaging techniques take two images andaverage to achieve better signal to noise ratio. In either of theseinstances the best image in terms of freedom from motion induced imageartefacts may be used. This could be chosen by a system operator oralternatively the best image could be chosen by suitable imageprocessing techniques. For instance an entropy based calculation couldbe conducted on the image to measure the amount of disorder of theimage, the image with the lower entropy having more order and lessartefacts. The skilled person would be well versed in suitable entropycalculation techniques.

[0086] Where no reference data is available a reference image could befabricated from the available data.

[0087]FIG. 4 shows a 2D slice of a subject. An image of the stationarysubject is shown on the left. This image was used as the referenceimage. A motion corrupted image was taken as shown in the centre (thesubject nodded during a scan). Using the method of the present inventiona corrected image was obtained as shown. It can be seen that thecorrected image is significantly improved over the motion corruptedimage.

[0088]FIG. 5 shows three orthogonal slices through a subject. Again areference image is shown on the left, a motion corrupted image shown inthe centre and an image corrected according to the present inventionshown on the right. Again an improvement can be seen in the correctedimage

[0089] In a further embodiment the cross-spectrum based method describedabove can be used to give an indication of confidence with which themotions have been determined. Under ideal conditions, where thereference data set and image data set differ only due to translationalmotion, the Fourier Transform of K_(c) is a peak with all other valueszero. In practice both data sets will be different due to noise andbecause they are acquired at different times and possibly underdifferent conditions. Therefore in practice the Fourier Transform willresult in a peak with other values that are non-zero. The ratio of theheight of the peak to the average of the other values can be used as ameasure of reliability, the higher the ratio the more reliable themotion determined. This could be used to simply give an indication ofconfidence in motion determination or could be used to indicate that adifferent reference data set should be tried. The method couldautomatically determine the ratio and produce a signal if the ratio fellbelow a certain threshold.

[0090] Other embodiments and advantages of the invention will beapparent to the skilled person without departing from the inventiondescribed herein.

1. A method of producing an image of a scanned object corrected forartefacts introduced by unwanted motion of said object during the scancomprising the steps of; taking a scan k-space data set derived from theobject scan, taking a reference k-space data set of the object,comparing the scan k-space data set and reference k-space data set todetermine the extent of any motion occurring during the scan in the scank-space data set, using the scan k-space data set and information on theextent of motion of the scan k-space data set to produce a motioncorrected scan k-space data set, and using the corrected scan k-spacedata set to produce an image.
 2. A method according to claim 1 whereinthe method includes the step of dividing the scan k-space data set andreference k-space data set into sections and wherein the step ofcomparing the scan k-space data set and reference k-space data setcomprises the step of comparing the scan k-space data set and referencek-space data set section by section to determine the extent of motion ofthe scan k-space data set.
 3. A method as claimed in claim 2 whereineach section of the scan k-space data set is a line of phase encodedata.
 4. A method as claimed in claim 2 wherein each section of the scank-space data set is a plane of k-space data.
 5. A method as claimed inany preceding claim wherein the method comprises a step of comparing thescan k-space data set and reference k-space data set to determine theextent of any rotation and separately comparing the scan k-space dataset and reference k-space data set to determine the extent of anytranslation.
 6. A method as claimed in any preceding claim wherein themethod includes the step of comparing the moduli of selected points inthe scan k-space data set and reference k-space data set to determinethe extent of rotational motion.
 7. A method according to claim 6wherein the selected points are the points with the largest modulusvalue in selected blocks of k-space.
 8. A method according to any ofclaims 6 or 7 wherein the step of comparing the moduli of selectedpoints comprises the steps of iteratively transforming the selectedpoints in either the reference k-space data set or scan k-space data setwith various possible rotational transformations and comparing with theother k-space data set.
 9. A method according to claim 8 wherein thetransformation is carried out using linear interpolation.
 10. A methodaccording to any preceding claim wherein a cross spectrum based methodis used to compare the scan k-space data set and reference k-space dataset to determine the extent of translational motion.
 11. A methodaccording to claim 10 wherein the cross-spectrum based method involvesthe step of combining the reference k-space data set and scan k-spacedata set to form a combined k-space data set, performing a Fouriertransform of the combined k-space data set and determining the positionand phase of the peak of the Fourier Transformed combined k-space dataset to determine the extent of translational motion.
 12. A methodaccording to claim 10 or claim 11 wherein in performing the crossspectrum analysis k-space is padded with complex values having a modulusequal to one and random phases.
 13. A method according to any of claims10-12 wherein the method includes the step of correcting the scank-space data set for the determined translational motion and comparingthe phase of the motion corrected scan k-space data set with the phaseof the reference k-space data set to determine sub-pixel motions.
 14. Amethod as claimed in any of claims 10-13 wherein in performing the crossspectrum analysis any data point known to be unreliable is replaced witha data point having a modulus equal to one and random phase.
 15. Amethod as claimed in claim 11 wherein, following the Fourier Transformof the combined k-space data set, the ratio of the peak amplitude toaverage amplitude of the rest of values is determined.
 16. A methodaccording to any preceding claim wherein the method includes an initialstep of aligning both the scan k-space data set and the referencek-space data set so that the maximum (DC) point coincides with thecentre of the data grid.
 17. A method according to any preceding claimwherein the method includes an initial step of globally scaling the twok-space data sets so at the maximum (DC) point the two k-space data setshave the same complex value.
 18. A method according to any precedingclaim wherein the images represented by the reference k-space data setand scan k-space data set are filtered to introduced a blur into thek-spaces before comparison of the k-space data sets.
 19. A methodaccording to any preceding claim wherein the scan k-space data set isobtained from a magnetic resonance imaging scanning means.
 20. A methodaccording to any preceding claim wherein parts of the image representedby the scan k-space data set is masked before comparison.
 21. An imagingsystem comprising a scanning means for scanning an object and deriving ascan k-space data set from signals received from the object and aprocessing means for manipulating the scan k-space data set to producean image of the scanned object corrected for artefacts introduced byunwanted motion of said object during the scan wherein the processingmeans is adapted to take the scan k-space data set derived from theobject scan, take a k-space reference data set of the object, comparethe scan k-space data set and reference k-space data set to determinethe extent of any motion in the scan k-space data set, use theinformation on the extent of motion of the scan k-space data set toproduce a motion corrected scan k-space data set, and use the motioncorrected scan k-space data set to produce an image.